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ABSTRACT 

Results and methods on three different areas of contemporary research are 
outlined. These include adaptive methods, the rolling contact problem for 
finite deformation of a hyperelastic or viscoelastic cylinder, and non- 
classical friction laws for modeling dynamic friction phenomena. 
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1 . INTRODUCTION 


This paper addresses three subjects that impact on the computer simula- 
tion of nonlinear tire behavior: adaptive methods, which represent schemes 

for assessing numerical error and automatically adapting the mesh so as to 
improve accuracy; the rolling contact problem, which is at the heart of tire 
analysis; and new friction laws, which are essential in developing realistic 
models of frictional contact. Space limitations preclude a detailed discus- 
sion of these issues; but further details can be found in recent papers and 
reports by the author and his colleagues [1-17]. 


2. ROLLING CONTACT 

The general rolling contact problem as a basis for nonlinear tire anal- 
ysis involves some of the most challenging and difficult problems in struc- 
tural mechanics. Among the complicating features are the presence of 
unilateral contact, friction, inertia effects, multi-parameter bifurcations, 
the emergence of standing waves, viscoelastic and thermal effects, large 
deformations, the necessity of modeling of complex materials such as fiber- 
reinforced rubbers, and the presence of non-conservative pressurization 
loadings. A first step toward resolving such issues is the formulation of 
correct kinematics and variational principles for a special case: the steady- 

state rolling of a hyperelastic or viscoelastic cylinder in contact with a 
rigid foundation and in a state of plane strain. 

The kinematical situation is shown in Fig. 1 where the geometry of the 
reference configuration (a rigid spinning cylinder with no contact) is com- 
pared to the geometry of a deformed cylinder in steady-state rolling contact 
with a rough (frictional) roadbed (foundation). 

Key features of the kinematic description are listed as follows: 

1) Time appears only implicitly in the formulation; if (R, (§), Z,) are 
material coordinates, the referential coordinates are 


r=R, 0 = @ + ut , z = Z 


where w is the angular velocity of the rigid, reference cylinder. 
2) If the motion is defined by the map 

x i = Xi (r, 0, z) i = 1, 2, 3 


where x^ = spatial Cartesian coordinates of particles in the current confi- 
guration, then velocity and acceleration components are 


3 X i 

v i = " JT 


a . “ 0) 
1 


2 
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3) The unilateral contact constraint can be expressed in the form 


X 2 < H on T c 


where is the exterior contact surface and (§) is the distance from the 

hub center to the foundation. This condition can also be written 


(X 2 ” H) + = 0 


where (•)+ denotes the positive part of (.). 

4) The time history of deformation can be expressed in terras of strains 
of particles located on the same circular arc in the reference configuration. 
For example, if E is the Green-Saint Venant strain tensor, its time history 
over an internal 0 < t < t satisfies: 

{E (r,0,z,x); 0 < t < t} = {e (r,x,z,t), 0 £ x £ 


This property makes it possible to incorporate viscoelastic effects into the 
rolling contact problem in a straightforward way. 

For illustration purposes, we consider a class of rolling contact prob- 
lems in which the following constitutive properties hold: 

a) The material is either an isotropic hyperelastic material char- 
acterized by a strain energy function 


W - WCIj.I^Ij) 

(or W = W(I^,I 2 ) if I 3 = 1 - an incompressible material) in which Ij,I 2 ,I 
are the principal invariants of the deformation tensor C = F T £, or the 
material is a viscoelastic material characterized by a linear viscoelastic 
perturbation of a hyperelastic material; e.g., 


S 


3W 

3E 


+ Y / 


k(T,t) E(x)dx 


with y a viscosity parameter, k a material kernel, and S the second 
Piola-Kirchhof f stress tensor. 



b) The normal stiffness of the contact interface obeys a constitutive 
law of the type 


m 

C n ' C „ (X 2 ' 10+ " 0 " r c 


where t n is the normal stress and c n and are material constants. 

If m n - 1 and c n = 1/e, where e is a positive constant, this relation 
coincides with the normal contact stress associated with an exterior penalty 
approximation of the unilateral constraint condition. 

c) If the cylinder is given a prescribed velocity v Q relative to 
the roadway, the slip velocity on the contact surface is 


W T = v o ” Xi = v Q - w9 0 x 


5) The friction law is (with the frictional stress) 


I~tI < = > w T = 0 

I t T ! = p|tj = > W T = -xj |t T | for some X > 0 


Variational principles for various rolling contact problems are summa- 
rized in Figs. 2-8, beginning with the pure spinning of a cylinder without 
contact and progressing to the general variational inequality for finite de- 
formation rolling contact with friction. Various spaces of admissible func- 
tions are defined in these figures as well as several nonlinear forms. In 
particular, 


A(x> n) 

A/ 

b(x> n) 

iV rsj 

c(x, n) 

/V (V 

Kx, n) 

/V A/ 

J(x, n) 
f(n) 


= the internal virtual work produced by the Piola-Ki rchhof f 
stresses T q 

= the virtual work produced by inertia (radial acceleration) 
effects per unit of angular velocity velocity ui 

= work terra due to normal compliance of the interface 

= a virtual work terra representing the work done by the 
hydrostatic pressure p (present when the material is 
incompressible) 

= the virtual work term due to frictional forces 

= the virtual work due to external forces 
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A finite-element code has been developed based on this general varia- 
tional principle which has the following features: 

1. Biquadratic (Q 2 ) elements are used to approximate the displacement 
field and, for incompressible materials, Pp discontinuous linear elements are 
used to approximate hydrostatic pressures 

2. The frictional functions are regularized in a standard way 

3. A Riks-Crisf ield method with Newton-Raphson corrections is used to 
solve the nonlinear systems of equations characterizing the discrete problem 

To date, an extensive set of numerical solutions has been obtained using 
these concepts and methods. Here only one example is cited, which is inter- 
esting because of the slow emergence of standing waves as the angular velocity 
is increased for a fixed penetration H of a hollow rubber cylinder into a 
rough roadway with coefficient of friction y = 0.03. Computed deformed 
shapes and stress contours are shown in Figs. 9-13 for various values of w. 

We notice the steady development of more-or-less periodic wavelets on the ex- 
terior surface which meet at points at which singularities appear to be 
forming. The presence of friction on the contact surface destroys the symme- 
try of this wavelet pattern. Mild viscoelastic effects, such as those in 
rubbers at moderate temperatures, do not appreciably alter the structure of 
these deformed geometries. 

The generality of the formulation and of the methods employed here makes 
it possible to study numerically a wide range of rolling contact problems. 
Further work shall involve generalizations of these results to three- 
dimensional rolling contact problems which include the effects of turning, 
steering forces, and tilting relative to the roadway plane. 


3. ADAPTIVE METHODS 

We shall now turn to the important subject of adaptive finite-element 
methods. Adaptive methods should have a significant impact on not only tire 
analysis but also on general computational structural mechanics in the rela- 
tively near future. 

In general, adaptive methods seek to change the structure of an approxi- 
mate method to improve the quality of the solution. By structure, we mean the 
mesh density, locations of nodes, and order of the local polynomials. By 
quality of an approximation, we generally mean the error in approximation in 
some appropriate norm. There are, thus, two primary aspects of any adaptive 
f i nit e-element method: 

1) The estimation of the error 

2) The adaptive strategy 

In the first of these, it is assumed that an initial approximation of the 
solution is known, perhaps from a computation on a coarse mesh, and that this 
rough solution can be used to obtain an a posteriori estimate of the local 
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error over each finite element. Once an estimate of the local error is known, 
one must call upon some technique to reduce the local error and thereby im- 
prove the quality of the solution. 

There are two general types of methods we have studied for a posteriori 
error estimation of the local error over each finite element. Once an esti- 
mate of the local error is known, one must call upon some technique to reduce 
the local error and thereby improve the quality of the solution. 

There are two general types of methods we have studied for a posteriori 
error estimation: 

1. Residual methods 

2. Interpolation (or a priori) methods 

As the name implies, residual methods make use of element residuals - the 
residual or unbalance left over when the approximate solution is substituted 
into the governing equations and edge conditions on each element. 

The residual itself (e.g., the equilibrium unbalance in element forces) 
is not necessarily a good indication of local error. Indeed, the local 
residual can be nearly zero while the error can be quite large. For this 
reason, it is generally necessary to calculate certain local error indicators 
<|> x which bound the error above and below in appropriate norms. The calcula- 
tion of error indicators generally requires the solution of special local 
problems over each element in which the element residuals enter as data. For 
example, in the model elliptic problem, 

-Au = f in ft CI|r2 


u = on 3(1 


(with A the Laplacian and f given), the finite-element solution u^ 
satisfies 


/ Vu. • Vv, dxdy = / fv,dxdy 
ft ft 

for arbitrary test function v^, and over each element K of a mesh, the 
residual is 

r h = - % - f 


274 



Over element K, an error indicator <j)^ is computed which satisfies 


/ • Vv dxdy 

K K 


9u h 

[ r, vdxdy + d> -r — vds 
K h 3K Sn 


for v in H^K). One can show that the error e, = u - u h in the energy norm 

21/2 n n 

( | I e h I I i q = {/ dxdy } ) satisfies the bound 



2 

i,n 


< I 

K 




2 

1,K 


Various residual methods differ in the way these error indicators are 
defined and calculated. There are some residual techniques which can produce 
sharp local error estimates in virtually any norm for certain classes of prob- 
lems. (See Demkowicz and Oden [4, 5]). These methods are not restricted to 
linear problems and have been used to produce error estimates in highly non- 
linear problems (see [7, 16]). 

The interpolation methods make use of the fact that the interpolation 
error over an element K of diameter h^ over an element on which poly- 
nomials of degree p are used is 

l U " n h U l 1,K = C \ P I u I 2 ,K 

where u is a given smooth function, H^u is its interpolant, 

Mi v = I |Vu| 2 dxdy 
’ K 

l u l2,K = { (u L * u xy + V )dxd >' 

and C is a constant independent of and u. The idea behind interpola- 

tion methods is to estimate | u [ 2 k using results of a coarse-mesh approxima- 
tion (e.g., using finite-difference methods or extraction methods [6]). The 
constant C cannot, in general, be determined, so such interpolation methods 
can only be used to assess relative error in various finite elements. 

Once an estimate of the error is obtained, the local error is reduced 
adaptively using one of the following techniques: 
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1. h-methods: the mesh size h is reduced and the number of elements 

is increased in regions of high error. 

2. p-methods: the mesh is fixed, but the local order p of the 

polynomial shape functions is increased. 

3. Moving raesh methods: the nodes in a finite-element mesh are moved 

and concentrated in areas of high error. 

4. Combined methods: these involve combinations of the above three 

techniques. 

We have developed test codes which employ all of these methods. The 
results of some tests are given in Figs. 14-20, and specific comments follow. 

1. In Fig. 14 we see a distorted mesh obtained using a moving mesh 
strategy on the driven cavity problem for an incompressible viscous fluid (see 
[7]). 

2. Figures 15, 16, and 17 contain computed results in adaptive schemes 
based on residual methods devised by Demkowicz and Oden [4, 5]. The results 
shown here are for transient heat conduction problems with dominate convection 
effects and for nonlinear Burgers 1 equation vector-valued solutions which 
simulate the Navier-St okes equations. A special h-method is used here which 
employs a fine grid and a coarse-grid approximation to estimate error. 

3. Results of a new p-method for Navier-Stokes equations in two- 
dimensions are illustrated in Figs. 18, 19 (see [1, 16]). Here a rather 
coarse mesh is used and errors are reduced by increasing local polynomial 
degrees from 1 to 2 to 3. Different shading in these figures indicates dif- 
ferent levels of local L^ error, with black cells indicating an error of less 
than 5 percent, grey an error of less than 10 percent, and white an error of 
over 20 percent. Such large local errors are reduced before the solution is 
allowed to advance in time. Remarkably, the so-called effectivity index 0 
for this problem (which represents the ratio of the estimated local error to 
the actual local error) for an L^ - norm varied from around 1.001 to 0.860 
for the time ranges considered in a test example. This suggests that 
residual-type error estimates based on p-type strategies can be very accurate, 
even for transient nonlinear problems on coarse meshes. 

Figure 20 shows refined raesh patterns for a class of linear elliptic 
problems in which a very fast vectorizable h-raethod is used in conjunction 
with an interpolation-type error estimator (see [6]). One interesting aspect 
of this work, indicated by different shadings of elements in the figure, is 
that the distinction between "optimal" meshes determined using a very sophis- 
ticated error estimator (see [17]) and very crude estimates ([6]) is negligi- 
ble whenever strong singularities are present. 


4. NON-CLASSICAL FRICTION LAWS 

In our most recent calculations of rolling contact, we have employed 
special interface constitutive equations for the normal compliance of the 
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interface and the tangential frictional forces. Some of these laws are men- 
tioned in Section 2 of this paper (see also Fig. 6). The physical interpreta- 
tion and the motivation of such models are discussed in references [14, 15, and 
18]. 
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v - SPACE OF ADMISSIBLE MOTIONS 


H ( n f V | n|( r 0 ) - Rj a.e.on r D } 

V ■ ( H = (n i/ite) | | J n( W(VT|) dv 0 1 < » ) 


CASE: W - POLYNOMIAL OF l c . Il c , lll c 

V - SOBOLEV SPACE OF ORDER ( 1 ,p) 

= (W 1 >P (Q 0 )) 2 Tor some p, 2 i p < oo 


CASE : FOR VISCOELASTIC MATERIAL 
U s V u [t\ € V| 

|J Qo F4 l r =-„(r.E(VTl)):VTl dv 0 |<oo) 


Figure 2 




VARIATIONAL PROBLEM .« BJIP 


VARIATIONAL P ROBLEM 

FIND A MOTION X c V SUCH THAT 
A(X,Tl) = <i> 2 b(X,T\) + f(T\) for all T\ € v 
BOUNDARY VALUE PROBLEM 

DivT 0 (X) + po &o = Po 9 e 2 X in Oo 

X " Ri on r D 

T 0 (X)no = 0 onr/r D = r c 

A(X.Tl) = J Qo OW(VX)/3VX) : Vt) dv 0 

or = 1 q 0 t oW : v- n <iv 0 
Tq(X) - VX S(X) = 1ST PIOLA-KIRCHHOFF STRESS 
B(X.'H) s /Q 0 3 e X 3 e ndv 0 

3 e x • 3 0 ti b ox/sexanj/ao) 

«'»i )= jQ 0 Po b o , n'iv 0 

Figure 3 

ROLLING CONTACT WITHOUT 
FRICTION 

% = UNILATERAL CONSTRAINT SET 
S ( t\ e v | t\ z <. H on r c ) 

VARIATIONAL INEQUALITY 

FIND A MOTION X e % SUCH THAT 

A(X.il-X) - w 2 B(X.il-X) * F(t|~X) for all T| e s( 

INCOMPRESSIBLE MATERIAL 


X= f v | n 2 i H on r c , detVTi = 1 a.e. in Q 0 } 

Figure 4 



ROLLING CONTACT WITH FRICTION 


<u = SPACE OF VELOCITY-MOTIONS 

= {t\€V \ (wOqTIi,^) £ v ) 

VIRTUAL POWER BY FRICTION 
j : v x a — > ^ 

j(X.Tl> s J rc JJ | T 0 CX)no n | |o)8 0 n,-v o | <Js 0 

VARIATIONAL INEQUALITY = B.V.P. 

VARIATIONAL INEQUALITY 

FIND A MOTION X f 11 n % SUCH THAT 

A(x.'n-v 0 , x) - u 2 B(x.‘n-v 0 , x) + j(x.'n) - 

j(X.V 9 , X)*f('n-V e , X) for all T| e j( 

with V 0 'X = (co9 0 X,,X2) 

Figure 5 


GENERALIZATION TO 
NONCLASS1CAL NORMAL 
CONTACT/FRICTION LAWS AND 
RFGlil ARIZFD FRICTION 
FUNCTIONALS 

NORMAL COMPLIANCE ON CONTACT SURFACE 
C : v » v — > Ji 

C(X.il) 2 J rc 3 e X2) ti 2 ds 0 
KX2.“3 e X2) = c n (X2’H) + m n + 
b n (X2 - H)+'n X 2 
X 2 = o)( 8X2/30) 

VARIATIONAL INEQUALITY WITH NORMAL COMPLIANCE 
FIND A MOTION X £ « SUCH THAT 

A(x.-n-v e 'X) - w^X.-n-Ve'x) + c(x,-n-v 0 'x) 
+ j(x.-n)-j(x.v 0 'x)2f(Ti-v 0 'x) 

Figure 6 


LAGRANGE MULTIPLIER METHOD 
3 = SPACE OF MULTIPLIER 
detV-q e LP(Q 0 ) — > Q = LP’(Q„) 

( 1 /p) + ( 1 +p) - 1 


FIND (X,p) 6 J (*« SUCH THAT 

A(x.ii-X) - w^cx.ti-x) + kp.x.ii-X) * ftii-x) 

for all T| e X 

( q,(detVX - 1 ) ) = 0 for all q e Q 


Kp.X.il) = TRILINEAR FORM 

= Jq o P ad JVX : V-q dv 0 
( q,(detvx- 1 ) ) = J Qo q (detVX - 1 ) dv 0 

Figure 7 


PFfil H API7FD FP1CT1QN LAW 

J f : IX * V -- > ^ 

J e (X.'n) s J rc p|T 0 (X)no n| 

V MeXi^o I )nids 0 
with <P f (| w t |) = 1 if |w t | > e 

= |w t |/e if|w t |se 

or = tanh( | w t | /() 

REMARK: lim f _ >0 J^X.H) - j(X.'H) “ KX.Vq'X) 

REGULARIZED VARIATIONAL EQUALITY 

FIND A MOTION X f f V , p € 9 SUCH THAT 

A(X e .H) - w^CXj.n) + CCXf.^l) + JfCXf.'H) + 
Kp.X e .^)- « ti) 

(q,det( 7 X f - 1))-0 


Figure 8 


for all n e v 
for all q t 9 


Cl =80 psi C2=20 psi VIS=0 FRIC=.3 DISP=R Q - H = .1 in. 



Figure 9 



Figure 10 
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Cl =80 psi C2 = 20 psi VIS=0 FRIC=.3 DISP=R 0 - H = .1 in. 



Figure 11 



.1 in. 


300 rpm 


Figure 12 
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Figure 13 



Driven cavity problem. Optimal mesh after 8 FE recalcula- 
tions. a = 4 , 6 = 4 , y = 0 


Figure 14 

C- 1 ^ 
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Heat equation with a dominating convection - 
solution after 1 time step, t =0.02 


Figure 15 




Heat equation with a dominating convection - 
solution after 25 time steps, t = 0.5 


Figure 16 



Burger's equations. First component of the 
solution after 1 time step, t= 0.02 

Figure 17 


♦ 




TIME=0. 1000E+00 



Problem with exponential solution. Deviations of the 

local effectivity ratios e? , from unity. Parameters: h = 2, 

At = 0.1, 6 = 0.25, M = 4 ^’ i 

Figure 18 
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TIME=0 . 5000E+00 



Problem with exponential solution. Deviations of the local 
effectivity ratios ejjj g from unity. Parameters h = 2, 

At = 0.1, 6 = 0.5, M =*4 

Figure 19 



1. Preconditioned Jacobi Conjugate Gradient Scheme 

2. 4-Elements/Refinement 

3. Rates = Uniform/Const. = (C) 1 ' 2 


Figure 20 


